If you're running an experiment on a population that is normally distributed with a known standard deviation and you observe a sample... what is the likelihood of that sample being the mean of the population you're observing?

That's what I'm solving for today.


In [40]:
import numpy, math

foo = [round(numpy.random.normal(20, 3),2) for x in range(0,5)]

In [41]:
foo


Out[41]:
[14.62, 24.55, 18.23, 21.83, 19.25]

How to calculate the distribution


In [42]:
n = len(foo)
sigma = 3
mu = 20

def likelihood_of_normal_parameters(data, mu, sigma):
    return numpy.exp(-n/(2*math.pow(sigma, 2))*math.pow(numpy.mean(data) - mu,2))

hypotheses = numpy.array([1]*10000)/(1.0*10000)

hypotheses


Out[42]:
array([ 0.0001,  0.0001,  0.0001, ...,  0.0001,  0.0001,  0.0001])

In [43]:
for i in range(0,10000):
    hypothesis = i/(1.0*100)
    hypotheses[i] = hypotheses[i]*likelihood_of_normal_parameters(foo, hypothesis, sigma)

In [44]:
new_hypotheses = hypotheses/numpy.sum(hypotheses)

In [45]:
import matplotlib.pyplot as plt

plt.plot(new_hypotheses)


Out[45]:
[<matplotlib.lines.Line2D at 0x10fe11fd0>]